// 输入文件：（当前目录下即可）
// line.dat: 要取索引点的位置坐标，　为线段，由pointwise输出，segment
// cavity_b.x: 网格坐标，为二进制Ｉ文件。
// 输出：
// index.dat: 位置坐标最接近的网格点索引
// testline.dat: 索引点对应的位置，用于确认是否为对应的点

//屏蔽visual studio vc 对标准c/c++ 语言的批判
#define _CRT_SECURE_NO_WARNINGS
#include <cmath>
#include <chrono>
#include <iostream>
#include <fstream>
// #include <vector>
// #include <valarray>
#include <string>
// #include <cstddef>
#include "multiarray.h"

// #define test
// #define HIGH
using namespace std;
struct Point
{
    double x, y, z;
};
struct Block
{
    int dimi, dimj, dimk;
    // double
    // Block()
};

inline istream &read_point(istream &, Point &);
inline ostream &print_point(ostream &, Point &);

int main()
{
    constexpr int three = 3;
    using namespace std::chrono;
    steady_clock::time_point t1 = steady_clock::now();

    // 读入线段，也就是感兴趣的点
    int npoint;
    string ifname1 = "line.dat";
    ifstream iline_file{ifname1, ios_base::in};

    iline_file >> npoint;
#ifdef HIGH
    MultiArray<Point> points(npoint);
    for (int i = 0; i != npoint; ++i)
    {
        iline_file >> points(i).x >> points(i).y >> points(i).z;
    }
#else
    const int xyz = three * npoint;
    double *points = new double[xyz];
    for (int i = 0; i != xyz; ++i)
    {
        iline_file >> points[i];
    }
#endif
    iline_file.close();

// 输出数据以便检查
#ifdef test
    for (int i = 0; i < npoint; ++i)
        print_point(cout, points(i));
    cout << endl;
#endif

//读入网格
#ifndef test
    string ifname2 = "cavity_b.x";
#else
    string ifname2 = "../../pointwise/cube/cube_b.x";
#endif // test
    ifstream mesh_file{ifname2, ios_base::in | ios_base::binary};
    int nblock;
    mesh_file.read(reinterpret_cast<char *>(&nblock), sizeof(int));
    cout << nblock << "\t blocks"
         << "\n";

    MultiArray<Block> blocks(nblock);
    mesh_file.read(reinterpret_cast<char *>(blocks.elem), sizeof(Block) * nblock);
#ifdef test
    for (int i = 0; i != nblock; ++i)
        cout << blocks(i).dimi << "\t" << blocks(i).dimj << "\t" << blocks(i).dimk << "\n";
        // cout << sizeof(Block) << "\n";
        // cout << sizeof(int) << "\n";
#endif // test

    // 每一块都是单独处理
    // 尽量使用内置类型，自己定义的再嵌套效率比较低，必要的使用自定义类型
    MultiArray<MultiArray<double>> grid(nblock);
    for (int i = 0; i != nblock; ++i)
    {
        grid(i).reconstruct(three, blocks(i).dimi, blocks(i).dimj, blocks(i).dimk);
        // cout << grid(i).size << "\n";
        // cout << grid(i).stride1 << "\n";
        // cout << grid(i).stride2 << "\n";
        // cout << grid(i).stride3 << "\n";

        mesh_file.read(reinterpret_cast<char *>(grid(i).elem), grid(i).size * sizeof(double));
        // cout << "success" << "\n";
    }
#ifdef test
    for (int i = 0; i != nblock; ++i)
    {
        for (int k = 0; k != grid(i).size; ++k)
            if (k % 5 != 0)
                cout << grid(i).elem[k] << "\t";
            else
                cout << "\n";
    }
#endif //test
    mesh_file.close();

    //距离最小
    double distance, tmp1;
    int tmp4, t5;
    MultiArray<int> index(2, npoint);
    for (int i = 0; i != npoint; ++i)
    {
        distance = 1000.0;
        t5 = i * 3;
        for (int j = 0; j != nblock; ++j)
        {
            tmp4 = grid(j).stride1;
            for (int k = 0, m = tmp4, n = tmp4 * 2; k != tmp4; ++k, ++m, ++n)
            {
                tmp1 = sqrt((grid(j).elem[k] - points[t5]) * (grid(j).elem[k] - points[t5]) +
                            (grid(j).elem[m] - points[t5 + 1]) * (grid(j).elem[m] - points[t5 + 1]) +
                            (grid(j).elem[n] - points[t5 + 2]) * (grid(j).elem[n] - points[t5 + 2]));
                if (tmp1 < distance)
                {
                    distance = tmp1;
                    index(0, i) = j;
                    index(1, i) = k;
                }
            }
        }
    }
    // cout << "hello" << endl;
#if 0
    for (int i = 0; i != npoint; ++i)
    {
        cout << index(0, i) +1 << "\t" << index(1, i) + 1 << "\n";
    }
#endif //test

    //转换与输出
    int ni, nj, nk;
    int sd1, sd2;
    ofstream fout{"multiarray_line_index2.dat", ios_base::out};
    fout << npoint << "\n";
    for (int i = 0; i != npoint; ++i)
    {
        sd1 = grid(index(0, i)).dim2 * grid(index(0, i)).dim3;
        sd2 = grid(index(0, i)).dim2;
        nk = index(1, i) / sd1;
        nj = (index(1, i) - nk * sd1) / sd2;
        ni = (index(1, i) - nk * sd1 - nj * sd2);
        fout << index(0, i) + 1 << "\t" << ni + 1 << "\t" << nj + 1 << "\t" << nk + 1 << "\n";
    }
    // cout << "ddddd" << "\n";
    steady_clock::time_point t2 = steady_clock::now();
    duration<double> dr = duration_cast<duration<double>>(t2 - t1);
    // printf("%lf\n", dr);;
    cout << dr.count() << endl;
    return 0;
}

istream &read_point(istream &is, Point &p)
{
    return is >> p.x >> p.y >> p.z;
}
ostream &print_point(ostream &os, Point &p)
{
    return os << p.x << "\t" << p.y << "\t" << p.z << "\n";
}